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Abstract 

We investigate the dynamics of matter-wave solitons in the presence of a spatially varying atomic 
scattering length and nonlinearity. The dynamics of bright and dark solitary waves is studied 
using the corresponding Gross-Pitaevskii equation. The numerical results are shown to be in very 
good agreement with the predictions of the effective equations of motion derived by adiabatic 
perturbation theory. The spatially dependent nonlinearity leads to a gravitational potential that 
allows to influence the motion of both fundamental as well as higher order solitons. 
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I. INTRODUCTION 



Recent years have seen enormous progress with respect to our understanding and the 
controlled processing of atomic Bose-Einstein condensates (BECs) [1] both for theory and 
experiment. In case of nonlinear excitations, specifically solitons, the experimental observa- 
tion of dark [2] , bright 0, 0] and gap [5] solitons has inspired many studies on matter- wave 
solitons in general. Apart from a fundamental interest in their behavior and properties, 
solitons are potential candidates for applications since there are possibilities to coherently 
manipulate them in matter- wave devices, such as atom chips j6]. Moreover, the formal sim- 
ilarities between matter-wave and optical solitons indicate that the former may be used in 
future applications similarly to their optical siblings, which have a time-honored history in 
optical fibers and waveguides (see, e.g., the recent reviews ll,!^). 

Typically dark (bright) matter- wave solitons are formed in atomic condensates with re- 
pulsive (attractive) interatomic interactions, i.e. for atomic species with positive (negative) 
scattering length a. One of the very interesting aspects for tailoring and designing the prop- 
erties of (atomic or molecular) BECs is the possibility to control the interaction of ground 
state species by changing the threshold collision dynamics and consequently changing either 
the sign or the magnitude of the scattering length. A prominent way to achieve this is to 
apply an external magnetic field which provides control over the scattering length because 
of the rapid variation in collision properties associated with a threshold scattering resonance 
being a Feshbach resonance (see refs. [HQ] and references therein). For low-dimensional 
setups a complementary way of tuning the scattering length or the nonlinear coupling at 
will is to change the transversal confinement in order to achieve an effective nonlinearity pa- 
rameter for the dynamics in e.g. the axial direction. In the limit of very strong transversal 
confinement this leads to the so-called confinement induced resonance at which the modified 
scattering length diverges [11]. A third alternative approach uses the possibility of tuning 
the scattering length with an optically induced Feshbach resonance [l2|. Varying the in- 
teractions and collisional properties of the atoms was crucial for a variety of experimental 
discoveries such as the formation of molecular BECs [13] or the revelation of the BEC-BCS 
crossover [14J. Recent theoretical studies have predicted that a time- dependent modulation 
of the scattering length can be used to prevent collapse in higher-dimensional attractive 
BECs [3], or to create robust matter- wave solitons flq . 
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Reflecting the increasing degree of control with respect to the processing of BECs it 
is nowadays "not only" possible to change the scattering length in the same way for the 
complete ultracold atomic ensemble i.e. it can be tuned globally, but it is possible to obtain 
a locally varying scattering length thereby providing a variation of the collisional dynamics 
across the condensate. According to the above this can be implemented by a (longitudinally) 
changing transversal confinement or an inhomogeneity of the external magnetic field in the 
vicinity of a Feshbach resonance. There exist only very few investigations on condensates in 
such an inhomogeneous environment ^ | . 

To substantiate the above, let us specify in some more detail the case of a magnetically 
tuned scattering length. The behavior of the scattering length near a Feshbach resonant 
magnetic field Bq is typically of the form a(B) = a[l — A/(B — Bq)], where a is the value of 
the scattering length far from resonance and A represents the width of the resonance (see, 
e.g., [3]). Let us consider a quasi-lD condensate along the x-direction exposed to a bias 
field B\ sufficiently far from the resonant value Bq in the presence of an additional gradient 
e of the field i.e. we have B = B\ + ex such that B\ > Bq + A (without loss of generality 
we take e > 0). Assuming exj[B\ — Bq) <C 1 for all values of x in the interval (— L/2, L/2), 
where L is the characteristic spatial scale on which the evolution of the condensate takes 
place, it is readily seen that the scattering length can be well approximated by the spatially 
dependent form a(x) = a + ctix, where a = a[l — A/(Bi — B )] and a x = eAa/(Bi — Bq) 2 . 
In the following we will assume that ao and a± are of the same sign. 

This opens the perspective of studying collisionally inhomogeneous condensates. In this 
work, we provide a first step in this direction by investigating the behavior of nonlinear 
excitations, specifically bright and dark matter-wave solitons in attractive and repulsive 
quasi one-dimensional (ID) BECs, in the presence of a spatially-dependent scattering length 
and nonlinearity. We investigate the soliton dynamics in different setups and analyze the 
impact of the spatially varying nonlinearity by numerically integrating the Gross-Pitaevskii 
GP) equation as well as in the framework of adiabatic perturbation theory for solitons 

mm 

The paper is organized as follows: In Sec. II the effective perturbed NLS equation is 
derived. In Sec. Ill, fundamental and higher-order soliton dynamics are considered and 
Bloch oscillations in the additional presence of an optical lattice are studied. Sec. IV is 
devoted to the study of dark matter-wave solitons, and in Sec. V the main findings of this 
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work are summarized. 



II. THE PERTURBED NLS EQUATION 

At sufficiently low temperatures, the dynamics of a quasi-one-dimensional BEC aligned 
along the x-axis, is described by an effective one- dimensional (ID) GP equation (see, e.g., 
2^ |) of the form: 

where ip(x,t) is the order parameter, m is the atomic mass, and V(x) is the exter- 
nal potential. Here we assume that the condensate is confined in a harmonic trap i.e., 
V(x) = {l/2)muj 2 x 2 where uj x is the confining frequency in the axial direction. The non- 
linearity coefficient g, accounting for the interatomic interactions, has an effective ID form, 
namely g = 2hau±, where u± is the transverse-confinement frequency and a is the atomic 
s-wave scattering length. The latter is positive (negative) for repulsive (attractive) conden- 
sates consisting of e.g. 87 Rb ( 7 Li) atoms. 

As discussed in the introduction we assume a collisionally inhomogeneous condensate i.e., 
a spatially varying scattering length according to a(x) = + a\X where a and a\ are both 
positive (negative) for repulsive (attractive) condensates. Moreover, if the characteristic 
length L for the evolution of the condensate implies \a%L\ < |ao|, it is readily seen that the 
function a(x) can be expressed as a(x) = sA(x), where A(x) = |ao| + \a>i\x is a positive 
definite function (for —L/2 < x < L/2) and s = sign(ao) = ±1 for repulsive and attractive 
condensates respectively. We can then reduce the original GP Eq. (JTJ to a dimensionless 
form as follows: x is scaled in units of the healing length £ = fij yJn g Q m, t in units of £/c 
(where c = ^/n g Q /m is the Bogoliubov speed of sound), the atomic density n = \ijj\ 2 is 
rescaled by the peak density no, and energy is measured in units of the chemical potential 
of the system « = gone,; in the above expressions go = 2haou± corresponds to the constant 
(dc) value ao of the scattering length. This way, the following normalized GP equation is 
readily obtained, 

ij£ = + V ^ + ^OOMV, (2) 

where V(x) = (l/2)Q 2 x 2 and the parameter Q = (2ao^o) ~ 1 (^x-/^±) determines the magnetic 
trap strength. Additionally, g(x) = 1 + Sx is a positive definite function and 6 = eA(Bi — 
Bo)- 1 [1 - A(Si - Bq)]' 1 is the gradient. 



Let us assume typical experimental parameters for a quasi- ID condensate containing iV ~ 
10 3 atoms and with a peak atomic density n$ « 10 8 m" 1 . Then, taking the scattering length 
a to be of the order of a nanometer we assume that the ratio of the confining frequencies 
lj x /ujj_ varies between 0.01 and 0.1. Therefore, the trap strength Q is typically O(10 -2 )- 
O(10 _1 ). Furthermore, we will assume that the field gradient e is also small and accounts 
for the leading order corrections of the gradient 8 in what follows. Thus, Q and 5 are the 
natural small parameters of the problem. 

We now introduce the transformation ip = uj ^/g to rewrite Eq. (0) in the following form: 

.du ld 2 u l2 7-,/ \ /o\ 

Apparently, Eq. (j3J) has the form of a perturbed NLS equation (of the focusing or defocusing 
type, for s = — 1 and s = +1 respectively), with the perturbation R(u) being given by 
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(4) 



The last two terms on the right hand side of Eq. Q are of higher order with respect 
to the perturbation parameter 5 than the second term and will henceforth be ignored (this 
will be discussed in more detail below). We therefore examine the soliton dynamics in the 
presence of the perturbation including the first two terms of Eq. (0J). 

III. BRIGHT MATTER- WAVE SOLITONS 
A. Fundamental solitons 

In the case s — — 1 and in the absence of the perturbation, Eq. Q represents the tra- 
ditional focusing NLS equation, which possesses a commonly known family of fundamental 



bright soliton solutions of the following form 



u(x, t) = r]sech[r](x — x )] exp[i(kx — 4>(t)} (5) 

where rj is the amplitude and inverse spatial width of the soliton, xo is the soliton center, 
the parameter k = dx^/dt defines both the soliton wavenumber and velocity, and finally 
4>{t) = (l/2)(fc 2 — r] 2 )t+(j) is the soliton phase (0 O being an arbitrary constant). Let us assume 



now that the soliton width rj~ l is much smaller than f2 _1//2 and 5~ l [namely the characteristic 
spatial scales of the trapping potential and the function g(x)] or, physically speaking, the 
potential V(x) and the function g(x) vary little on the soliton scale. In this case, we may 
employ the adiabatic perturbation theory for solitons j^cj to treat analytically the effect of 
the perturbation R(u) on the soliton (J3J. According to this approach, the soliton parameters 
7], k and Xq become unknown, slowly-varying functions of time t, but the functional form of 
the soliton (see Eq. ©) remains unchanged. Then, from Eq. (J3J), it is found that the number 
of atoms N = f 00 \u\ 2 dx and the momentum P = (i/2) f_ °° [u(du* / dx) — u*(du/dx)] dx 
which are integrals of motion of the unperturbed system, evolve, in the presence of the 
perturbation, according to the following equations, 



dN 

— = -21m 
dt 

dP 

— = 2Re 

dt 



+00 

Ru*dx 



DC 
OO 



OO 



R——dx 

ox 



(6) 
(7) 



We remark that the number of atoms iV is conserved for Eq. (J2J) but the transformation, 
leading to Eq. Q no longer preserves that conservation law, leading, in turn, to Eq. (JBJ). 

We now substitute the ansatz (jSJ) (but with the soliton parameters being functions of 
time) into Eqs. (|B|) -([?J): furthermore we use a Taylor expansion of the second term of Eq. 

around x = xq (keeping the two leading terms). The latter expansion is warranted by 
the exponentional localization of the wave around x = xq. We then obtain the evolution 
equations for rj (t) and k(t), 

(h ' k V ^-ln(g), (8) 



dt dx 
dk dV rf d 



Ma)- (9) 



dt dxo 3 dxo 

To this end, recalling that dxo/dt = k, we may combine Eqs. ©-© to derive the following 
equation of motion for the soliton center: 

rf 2 x _ dV | T] 2 (0) fdg*\ (lQ) 



dt 2 dx Qg 2 (0) \dxo / 
where 77(0) and g(0) = g(x (0)) are the initial values of the amplitude and function g(x) 
respectively. Notice that the above result indicates that the main contribution from the 
spatially dependent scattering length comes to order 5 2 (while the contribution of the last 
two terms in Eq. (jlj) would have been 0(<5 3 ) and is neglected). It is clear that in the 



particular case where g(x) = 1 + 8x, Eq. (JTUj) describes the motion of a unit mass particle 
in the presence of the effective potential 

V eS (x ) = -Lol s x 2 Q - /3x Q , (11) 

where the parameter (3 is defined as 

0. „ (12) 

3[l + fa„(0)] 2 ' 

and 

^ bs = v 7 ^ 2 - W, (13) 

is the oscillation frequency of the bright soliton. In the absence of the spatial variation of 
the scattering length (5 = 0), Eq. (jlUJl actually expresses the Ehrenfest theorem, implying 
that the bright soliton oscillates with a frequency = ^ in the presence of the harmonic 
potential with strength Q. Nevertheless, the presence of the gradient modifies significantly 
the bright soliton dynamics as follows: First, as seen by the second term in the right-hand 
side of Eq. (J11J) . apart from the harmonic trapping potential, an effective gravitational 
potential is also present, which induces an acceleration of the initial soliton towards larger 
values of Xq (for 5 > 0) i.e. it shifts the center of the harmonic potential from xq = 
to Xo = At- Second, the oscillation frequency of the bright soliton is modified for 5^0, 

^bs 

according to Eq. (|13p. Moreover, depending on the initial values of the parameters, namely 
for (35 > Q 2 an interesting situation may occur, in which the effective harmonic potential, 
instead of being purely attractive, it can effectively disappear, or be expulsive. 

The solution to Eq. (|1(J|) in the variables y = x — (3 /u 2 s is, of course, a simple classical 
oscillator 

Vo{t) = 2/o(0) cos(u bs t) + ^JL S m(uj bs t) (14) 

^bs 

which is valid for u 2 s > 0. For uj 2 s < the trigonometric functions have to be replaced by 
hyperbolic ones. In the case u 2 s = 0, the resulting motion (to the order examined) is the 
one due to a uniform acceleration with x (t) = x (0) + x (0)t + (3t 2 /2. 

The above analytical predictions have been confirmed by direct numerical simulations. 
In particular, we have systematically compared the results obtained from Eq. (|1(J|) with 
the results of the direct numerical integration of the GP Eq. @. In the following, we use 
the trap strength Q = 0.05, initial soliton amplitude r/(0) = 1 and initial location of the 




FIG. 1: Top panels: The effective potential V(xo) as a function of the soliton center xq for a trap 
strength Vt = 0.05 for a fundamental bright soliton of amplitude 7/(0) = 1, initially placed at the 
trap center (xo(0) = 0). Different values of the gradient modify the character of the potential: 
in panel (a) it is purely attractive (S = 0, Q/2, \/2tl), while in (b) it is either purely gravitational 
(e) = \/3^) or expulsive (5 = 2S7). Bottom panels: Evolution of the center of the bright soliton for 
the above cases: (c) for attractive effective potentials and (d) for gravitational or expulsive ones. 
The agreement between numerical results (solid and dashed lines) and the theoretical predictions 
(triangles, dots) is excellent. 

soliton xo(0) = 0, and different values for the normalized gradient 8. The above values of the 
parameters, may correspond to a 7 Li condensate containing iV « 4000 atoms, confined in a 
quasi- ID trap with frequencies uo x = 2ir x 14 Hz and u± = 100u x . Note that these values 
correspond to a scattering length a = —0.21 nm (pertaining to a magnetic field 425 Gauss), 
a value for which a bright matter- wave soliton has been observed experimentally (^J . 

In Fig. [Ha), the original harmonic trapping potential V(x) (dotted line, 5 = 0) is 
compared to the effective potential modified by the presence of the gradient for 5 = (1/2)Q 
(dashed line) and 5 = \/2VL (solid line). In these cases, the effective harmonic potential is 
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attractive and the gradient displaces the equilibrium point to the right. Fig. IHb) shows the 
case 5 = \^3£l (dashed line) for which the effective harmonic potential is canceled resulting 
in a purely gravitational potential. Also, upon suitably choosing the value of 5, e.g., for 
5 = 2Q (solid line) the effective potential becomes expulsive. The dynamics of the bright 
matter-wave soliton pertaining to the above cases are shown in Figs. QJc) (for attractive 
effective potential) and^d) (for gravitational or expulsive effective potential). In particular, 
in Fig. [He), it is clearly seen that the evolution of the soliton center xq is periodic, but 
with a larger amplitude and smaller frequency of oscillations, as compared to the respective 
case with 5 = 0. The analytical predictions of Eq. (jl()j) - (jl4jl (triangles for 5 = y/20, and 
dots for 5 = (1/2)Q) are in perfect agreement with the respective results obtained by direct 
numerical integration of the GP Eq. On the other hand, as shown in Fig. CJd), in the 
case of a gravitational or expulsive effective potential, the function xo(t) is monotonically 
increasing, with the analytical predictions being in excellent agreement with the numerical 
simulations. For a purely gravitational or expulsive effective potential, Eq. (|SJ) shows that 
the amplitude (width) of the soliton increases (decreases) monotonically as well, which 
recovers the predictions of Ref. [l2|- This type of evolution suggests that the bright soliton 
is compressed adiabatically in the presence of the gradient. 

Let us consider another setup which combines the "effective" linear potential with an 
external harmonic and a periodic trap: 

V(x) = -fiV + V sin 2 (/tx) (15) 

The periodic potential in Eq. (fTo^) can be obtained experimentally by superimposing two 
counter-propagating laser beams. It is well-known that the dynamics in the combined pres- 
ence of a(n effective) linear and a periodic potential results in the so-called Bloch oscillations 
(for a recent discussion of the relevant phenomenology and bibliography see e.g. 22]). These 
oscillations occur due to interplay of the linear and periodic potential with a definite period 
T = 2k/ 13 22j. We have examined numerically this analytical prediction in the presence of 



an optical lattice potential with Vq = 0.25 and k = 0.5. The numerical evaluation of the 
period of the soliton motion in the combined potential is T « 22.15 less than 4% off the 
corresponding theoretical prediction. The time-periodic evolution of the soliton is shown in 
the spatio-temporal contour plot of Fig. |21 
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FIG. 2: Spatio-temporal contour plot (of the wavefunction square modulus) of a solitary wave for 
O = 0.075, 5 = \/30 and an optical lattice with V$ = 0.25 and k = 0.5. One can clearly discern 
the presence of Bloch oscillations in the evolution of the density, whose period is in very good 
agreement with the corresponding theoretical prediction. 



B. Higher-order solitons 



Apart from the fundamental bright soliton in Eq. (jHJ), it is well known 25j that specific 
AA-soliton exact solutions in the unperturbed NLS equation (Eq. (jSJ) with s = — 1 and 
R = 0) are generated by the initial condition u(x, 0) = v4sech(x) (for 77 = 1), and the soliton 
amplitude A is such that A — 1/2 < 7V< A + 1/2 to excite a soliton of order M . The exact 
form of the A/"-soliton is cumbersome and will not be provided here; nevertheless, it is worth 
noticing some features of these solutions: First, the number of atoms of the jV-th soliton 
is Af 2 larger than the one of the fundamental soliton and second, for any TV, the soliton 
solution is periodic with the intrinsic frequency of the shape oscillations being co> intr = 4rj 2 . 
We now examine the dynamics of the jV-soliton solution in the presence of the spatially 
varying nonlinearity. 

We have performed numerical simulations in the case of the so-called double (J\f— 2) 
bright soliton solution with initial soliton amplitude A = 2.5. In the absence of the gradient 
(5 = 0), if the soliton is placed at the trap center (xq = 0, with xq being the soliton 
center), it only executes its intrinsic oscillations with the above mentioned frequency uj intT . 
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On the other hand, if the soliton is displaced (xq ^ 0), apart from its internal vibrations, 
it performs oscillations governed by the simple equation oto + Q 2 xq = 0, in accordance to 



the Kohn theorem (see 261 and 



271 ] for an application in the context of bright matter- wave 



solitons). Nevertheless, for 5^0, the double soliton (initially placed at the trap center), 
contrary to the previous case, splits into two single solitons, with different amplitudes due 
to the effective gravity discussed in the case of the fundamental soliton. Due to the effective 
gravitational force, the soliton moving to the right (see, e.g., Fig. |3J) is the one with the 
larger amplitude (and velocity) and is more mobile than the one moving to the left (which 
has the smaller amplitude). 

As each of these two solitons is close to a fundamental one, their subsequent dynamics 
(after splitting) may be understood by means of the effective equations of motion derived 
in the previous section. In particular, depending on the values of the relevant parameters 
involved in Eq. [i](0) is now the amplitude of each soliton after splitting] the solitons 
may both be trapped, or may escape (either one or both of them), if the effective potential 
is expulsive. In the former case, both solitons perform oscillations (in the presence of the 
effective attractive potential) and an example is shown in Fig. El (for Q = 0.1, 5 = 0.01). 
Note that the center of mass of the ensemble oscillates with a period T = 2tt/Q = 62.8 
(which is in accordance with Kohn's theorem). During the evolution, as each of the two 
solitons oscillate in the trap with different frequencies, they may undergo a head-on collision 
(see, e.g., bottom panel of Fig. El at t « 55). It is clear that such a collision is nearly elastic, 
with the interaction between the two solitons being repulsive. 

Importantly, for smaller values of the trap strength Q, we have found that it is possible to 
release either one or both solitons from the trap: In particular, for Q = 0.05 and 5 = 0.025, 
we have found that the large amplitude soliton escapes the trap, while the small-amplitude 
one performs oscillations. On the other hand, for the same value of the trap strength but 
for 5 = 0.05 both solitons experience an expulsive effective potential and thus both escape 
from the trap. It is therefore in principle possible to use the spatially varying nonlinearity 
not only to split a higher-order bright soliton to a chain of fundamental ones, but also to 
control the trapping or escape of the resulting individual fundamental solitons. 
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FIG. 3: Evolution of a double soliton with initial amplitude A = 2.5 initially placed at the trap 
center (x = 0) with a strength $7 = 0.1 in the presence of a gradient 5 = 0.01. Top panel: Spatio- 
temporal contour plot of the density. Bottom panels: Snapshots of the evolution of the density are 
shown for t = (top left panel), t = 10 (top right panel), t = 30 (bottom left panel), and t = 60 
(bottom right panel), covering almost one period of the oscillation. Dashed lines correspond to the 
trapping potential. 

IV. DARK MATTER- WAVE SOLITONS 



We now turn to the dynamics of dark matter-wave solitons in the framework of Eq. (J2J 
for s = +1 (i.e., the defocusing case of condensates with repulsive interactions). Firstly we 
examine the equation governing the background wavefunction. The latter is taken in the 
form = $(x) exp(-ifit) (/i being the chemical potential) and the unknown background 
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wave function satisfies the following real equation, 

^ + \^-g{x)& = V{x)<$>. (16) 

To describe the dynamics of a dark soliton on top of the inhomogeneous background satis- 
fying Eq. (fTHjl. we introduce the ansatz (see, e.g., j^) 

ip — $(x) exp(— i/it)v(x, t), (17) 

into Eq. (j2J, where the unknown wavefunction v(x,t) represents a dark soliton. This way, 
employing Eq. (fTEj). the following evolution equation for the dark soliton wave function is 
readily obtained: 

Taking into account that in the framework of the Thomas-Fermi approximation jj| a simple 
solution of Eq. (fTtjj) is expressed as 



<I>U) |/max{ M ^ (X) ,0}, (19) 
equation (fTHj) can be simplified to the following defocusing perturbed NLS equation, 

^§-*f-i)-eH (20) 

where the perturbation Q(v) has the form, 

Q( v ) = (1 - H 2 ) + 



2(/x — V) dx dx 

and higher order perturbation terms have once again been neglected. In the absence of the 
perturbation, Eq. (|2*U|l represents the completely integrable defocusing NLS equation, which 



has a dark soliton solution of the form 



29| (for/i = l), 



v(x, t) = cosy?tanh£ + i simp, (22) 

where ( = cos<^ [x — (sin </>)£], while cosy? and sinyj are the soliton amplitude and velocity 
respectively, ip being the so-called soliton phase angle (\(f)\ < vr/2). To treat analytically the 
effect of the perturbation (|2*T|) on the dark soliton, we employ the adiabatic perturbation 
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theory devised in Ref. 2l|. As in the case of bright solitons, according to this approach, the 
dark soliton parameters become slowly-varying unknown functions of t, but the functional 
form of the soliton remains unchanged. Thus, the soliton phase angle becomes ip — > ip(t) 
and, as a result, the soliton coordinate becomes ( — > ( = cos ip(t) [x — x (t)], where 

x (t) = [ sin (p(t')dt', (23) 



is the soliton center. Then, the evolution of the parameter p> governed by the equation 



2l|, 



dip 
~dt 



1 



-Re 



+00 dv * ■ 
Q(v)—dx 



(24) 



2 cos 2 ip sin ip 

leads (through similar calculations and Taylor expansions as for the bright case) to the 
following result: 



— = — COS (D 

dt r 



ldV 1 d , , ,' 
+ -—Hg) 



2dx Q ' Sdx ' ' 

To this end, combining Eqs. (|2*H|) and (|23|) . we obtain the corresponding equation of motion 
for the soliton center, 

d 2 x ldV Id 



(26) 



dt 2 2 dx 3 dx 

in which we have additionally assumed nearly stationary dark solitons with cos<£> ~ 1. As 
in the case of bright solitons, the validity of Eq. (|26|) does not rely on the specific form of 
g(x), as long as this function (and the trapping potential) are slowly- varying on the dark 
soliton scale (i.e., the healing length). In the particular case with g(x) = 1 + 8x, Eq. (|2T)j) 
describes the motion of a unit mass particle in the presence of the effective potential 



WcfrOro) = ^ 2 x 2 + - ln(l + 6x ). 



(27) 



For 5 = Eq. (j2*U|) implies that the dark soliton oscillates with a frequency Vt/\/2 in the 
harmonic potential with strength Q 2^13- However, in the presence of the gradient, and 
for sufficiently small 5, Eq. ()27|) implies the following: First, the oscillation frequency of 
the dark soliton is downshifted in the presence of the linear spatial variation of the scattering 
length, according to 

(28) 



-Q 2 - -5 2 . 
2 3 



Additionally to the effective harmonic potential, the dark soliton dynamics is also modified 
by an effective gravitational potential (~ 5x /3), which induces an acceleration of the soliton 
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FIG. 4: (a) Two snapshots of the density of the dark soliton (at t = and t = 90) on top 
of a Thomas-Fermi cloud. The chemical potential is fi = 1, the trap strength is $7 = 0.05 and 
the gradient is 6 = 0.01. (b): The motion of the center of a dark soliton. Solid line and dots 
respectively correspond to the numerical integration of the GP equation and analytical predictions 
[see Eq. (|26|)]. respectively. 

towards larger values of Xq (for 5 > 0). It should be noted that as dark solitons behave as 
effective particles with negative mass, the effective gravitational force possesses a positive 
sign, while in the case of bright solitons (which have positive effective mass) it has the usual 
negative sign [see Eqs. (|TT|) and (|27j)]. 

Direct numerical simulations confirm the above analytical findings. In particular, we 
consider an initially stationary dark soliton (with cos<£>(0) = 0), placed at xq = 0, on top 
of a Thomas-Fermi cloud [see Eq. (|19|)] characterized by a chemical potential \x = 1 (the 
trapping frequency is here Q = 0.05). In the absence of the gradient such an initial dark 
soliton should be purely stationary. However, considering a gradient with 5 = 0.01, it is clear 
that the TF cloud will become asymmetric, as shown in Fig. Ufa) and the soliton will start 
performing oscillations. The latter are shown in Fig. HJ where the analytical predictions 
(points) are directly compared to the results obtained by direct numerical integration of 
the GP equation (solid line). As it is seen, the agreement between the two is very good; 
additionally, we note that the oscillation frequency found numerically is 27r/177, while the 
respective theoretical prediction is 27r/180.7, with the error being m 3%. 
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V. SUMMARY 



We have analyzed the dynamics of bright and dark matter- wave solitons in quasi- ID 
BECs characterized by a spatially varying nonlinearity. The formulation of the problem is 
based on a Gross-Pitaevskii equation with a spatially dependent scattering length induced 
e.g. by a bias magnetic field near a Feshbach resonance augmented by a field gradient. The 
GP equation has been reduced to a perturbed nonlinear Schrodinger equation, which is then 
analyzed in the framework of the adiabatic approximation in the perturbation theory for 
solitons, treating them as quasi-particles. This way, effective equations of motion for the 
soliton centers (together with evolution equations for their other characteristics) were derived 
analytically. The analytical results were corroborated by direct numerical simulations of the 
underlying GP equations. 

In the case of bright matter-wave solitons initially confined in a parabolic trapping po- 
tential, it is found that (depending on the values of the gradient and the initial soliton 
parameters), there is a possibility to switch the character of the effective potential from 
attractive to purely gravitational or expulsive. It has been thus demonstrated that a bright 
soliton can escape the trap and be adiabatically compressed. On the other hand, considering 
the additional presence of an optical lattice potential, it has been shown that in the case 
where the effective potential is purely gravitational, Bloch oscillations of the bright solitons 
are possible. Higher-order bright solitons have been shown to typically split in the pres- 
ence of a spatially varying nonlinearity to fundamental ones, whose subsequent dynamics 
is determined by the properties of the resulting single-soliton splinters. In the case of dark 
matter- wave solitons, the relevant background, i.e., the Thomas- Fermi cloud is modified by 
the inhomogeneous nonlinearity. The dynamics of the dark solitons follows a Newtonian 
equation of motion for a particle with a negative effective mass and the oscillation frequency 
of the dark solitons has been derived analytically. The latter is always down-shifted as com- 
pared to the oscillation frequency pertaining to a spatially constant scattering length. Thus, 
generally speaking, the presented results show that a spatial inhomogenity of the scattering 
length induced e.g. by properly chosen external magnetic fields is an effective way to control 
the dynamics of matter-wave solitons. 
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